#include <kuai/sim/nonbond.hpp>
#include <kuai/sim/efunctor.hpp>

#define TRY_TO_CREATE(function)												\
	if (name == boost::algorithm::to_lower_copy(String( #function ))) {		\
		if (hasCharge) {													\
			if (pbc) {														\
				return new NonbondEnergy<Functor##function, true, true>();	\
			}																\
			else {															\
				return new NonbondEnergy<Functor##function, true, false>();	\
			}																\
		}																	\
		else {																\
			if (pbc) {														\
				return new NonbondEnergy<Functor##function, false, true>();	\
			}																\
			else {															\
				return new NonbondEnergy<Functor##function, false, false>();\
			}																\
		}																	\
	}																		\

namespace kuai {

	BasicNonbondEnergy* createNonbondEnergyFunction(const char szName[], const PBC* pbc, bool hasCharge)
	{
		String name = splitPair(szName, '/').second;
		boost::algorithm::to_lower(name);

		TRY_TO_CREATE(LJ12_6);
		TRY_TO_CREATE(LJ_9_6);
		return NULL;
	}

	namespace {
	
		void addLJ12_6(RealNumber r2, const RealNumber par[], RealNumber e[]) {
			RealNumber rr2 = par[0] * par[0] / r2;
			RealNumber rr6 = rr2 * rr2 * rr2;
			e[0] += par[1] * rr6 * (rr6 - 2);
		}

		void addLJ12_6(XYZ& delta, RealNumber r2, Index i1, Index i2, RealNumber par[], RealNumber e[], XYZ gred[]) {
			RealNumber rr2 = 1 / r2;
			RealNumber rr6 = par[0] * par[0] * r2;
			rr6 = rr6 * rr6 * rr6;
			e[0] += par[1] * rr6 * (rr6 - 2);

			delta *= 12 * par[1] * rr6 * (1 - rr6) * rr2;
			gred[i1] += delta;
			gred[i2] -= delta;
		}
	}

	template <>
		void NonbondEnergy<FunctorLJ12_6, false, true>::addEItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[]) const
	{
		addLJ12_6(dot(delta, delta), pars+iPar, e);
	}
	template <>
		void NonbondEnergy<FunctorLJ12_6, false, false>::addEItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[]) const
	{ 
		addLJ12_6(dot(delta, delta), pars+iPar, e);
	}

	template <>
		void NonbondEnergy<FunctorLJ12_6, false, true>::addEGItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[], XYZ gred[]) const
	{
		addLJ12_6(delta, dot(delta, delta), i1, i2, pars+iPar, e, gred);
	}
	template <>
		void NonbondEnergy<FunctorLJ12_6, false, false>::addEGItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[], XYZ gred[]) const
	{
		addLJ12_6(delta, dot(delta, delta), i1, i2, pars+iPar, e, gred);
	}

	template <>
		void NonbondEnergy<FunctorLJ12_6, false, true>::addE_ABC(SimulationEnvironment& env, RealNumber e[])
	{
		if (needUpdateABCBuffer(env)) {
			updateABCBuffer(env);
		}
		std::pair<size_t, size_t> range = env.range(bufferList.size()/3);
		RealNumber cutoff2=cutoff*cutoff;
		for (size_t i = range.first; i < range.second; ++i) {
			Index atom1 = bufferList[i*3];
			Index atom2 = bufferList[i*3+1];
			Index iPar  = bufferList[i*3+2];
			XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
			pbc->update(delta);
			RealNumber r2;
			if (shorter2(delta, cutoff2, r2)) {
				addLJ12_6(r2, pars+iPar, e);
			}
		}
	}
	template <>
	void NonbondEnergy<FunctorLJ12_6, false, false>::addE_ABC(SimulationEnvironment& env, RealNumber e[]) {
		if (needUpdateABCBuffer(env)) {
			updateABCBuffer(env);
		}
		std::pair<size_t, size_t> range = env.range(bufferList.size()/3);
		RealNumber cutoff2=cutoff*cutoff;
		for (size_t i = range.first; i < range.second; ++i) {
			Index atom1 = bufferList[i*3];
			Index atom2 = bufferList[i*3+1];
			Index iPar  = bufferList[i*3+2];
			XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
			RealNumber r2;
			if (shorter2(delta, cutoff2, r2)) {
				addLJ12_6(r2, pars+iPar, e);
			}
		}
	}

	template <>
		void NonbondEnergy<FunctorLJ12_6, false, true>::addEG_ABC(SimulationEnvironment& env, RealNumber e[], XYZ gred[])
	{
		if (needUpdateABCBuffer(env)) {
			updateABCBuffer(env);
		}
		std::pair<size_t, size_t> range = env.range(bufferList.size()/3);
		RealNumber cutoff2=cutoff*cutoff;
		for (size_t i = range.first; i < range.second; ++i) {
			Index atom1 = bufferList[i*3];
			Index atom2 = bufferList[i*3+1];
			Index iPar  = bufferList[i*3+2];
			XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
			pbc->update(delta);
			RealNumber r2;
			if (shorter2(delta, cutoff2, r2)) {
				addLJ12_6(delta, r2, atom1, atom2, pars+iPar, e, gred);
			}
		}
	}
	template <>
		void NonbondEnergy<FunctorLJ12_6, false, false>::addEG_ABC(SimulationEnvironment& env, RealNumber e[], XYZ gred[])
	{
		if (needUpdateABCBuffer(env)) {
			updateABCBuffer(env);
		}
		std::pair<size_t, size_t> range = env.range(bufferList.size()/3);
		RealNumber cutoff2=cutoff*cutoff;
		for (size_t i = range.first; i < range.second; ++i) {
			Index atom1 = bufferList[i*3];
			Index atom2 = bufferList[i*3+1];
			Index iPar  = bufferList[i*3+2];
			XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
			RealNumber r2;
			if (shorter2(delta, cutoff2, r2)) {
				addLJ12_6(delta, r2, atom1, atom2, pars+iPar, e, gred);
			}
		}
	}

}
